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Abstract 

A Fortran 90 module GammaCHI for computing and inverting the gamma 
and chi-square cumulative distribution functions (central and noncentral) is 
presented. The main novelty of this package are the reliable and accurate 
inversion routines for the noncentral cumulative distribution functions. Ad¬ 
ditionally, the package also provides routines for computing the gamma func¬ 
tion, the error function and other functions related to the gamma function. 
The module includes the routines cdfgamC, invcdfgamC, cdfgamNC, 
invcdfgamNC, errorfunction, inverfc, gamma, loggam, gamstar and 
quotgamm for the computation of the central gamma distribution function 
(and its complementary function), the inversion of the central gamma dis¬ 
tribution function, the computation of the noncentral gamma distribution 
function (and its complementary function), the inversion of the noncentral 
gamma distribution function, the computation of the error function and its 
complementary function, the inversion of the complementary error function, 
the computation of: the gamma function, the logarithm of the gamma func¬ 
tion, the regulated gamma function and the ratio of two gamma functions, 
respectively. 
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small input function values P^(x,y), Qn(x,y) (lower than 10 -150 ) are not admis¬ 
sible. 

The admissible input parameter ranges for computing the noncentral cumu¬ 
lative gamma distribution functions P^(x,y), Q^(x,y) in standard IEEE double 
precision arithmetic are 0 < x < 10000, 0 < y < 10000, 0.5 < y < 10000 and 
the related parameter ranges for the noncentral chi-square cumulative distribution 
function. In the inversion of the noncentral gamma/chi-square distribution func¬ 
tions, very small input function values P^(x,y)^ Q^(x,y) (lower than 10 -25 and 
Iq- 35, respectively) are not admissible. 

Running time: 

It varies depending on the function and the parameter range. 


1. Introduction 

The Fortran 90 module GannnaCHI provides reliable and fast routines 
for the inversion and computation of the gamma and chi-square distribution 
functions. These functions appear in many problems of applied probability 
including, of course, a large number of problems in Physics. The module 
also includes routines for the computation of the gamma function, the error 
function and its complementary function, the inverse of the complementary 
error function, the logarithm of the gamma function, the regulated gamma 
function and the ratio of two gamma functions. 

The main novelty of the package GammaCHI are the inversion routines 
for the noncentral cumulative distribution functions. These algorithms are 
based on the asymptotic methods presented in [lj although improvements in 
the estimation of the initial values for small values of the parameter /i are 
included in the present algorithms. Among other applications, the compu¬ 
tation and the inversion of the noncentral cumulative distribution functions 
included in GammaCHI appear in the analysis of signal detection in differ¬ 
ent physical scenarios such as optics (2|, radiometry [3| or quantum detection 
Q, |f|. Also, the inversion routines can be used to generate normal variates 
as well as central and noncentral gamma (or chi-square) variates, which can 
be useful, for example, in Monte Carlo simulations. 

2. Theoretical background 

2.1. Central gamma distribution function 

The incomplete gamma functions are defined by 
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•X 


•oo 


( 1 ) 


7 (a, x) = / t a 1 e 1 dt, T(a, x) = / t a 1 e t dt. 


Let us define 


P{a , x ) = x), Q(a, x ) = (a, x), 

i (a) i (a) 


( 2 ) 


where we assume that a and x are real positive numbers. 

The functions P(a,x), Q(a,x) are the central gamma distribution func¬ 
tion and its complementary function, respectively. These distribution func¬ 
tions which appear in many problems of applied probability include, as partic¬ 
ular cases, the standard chi-square probability functions P(x 2 \ v ) and Qix 2 ^) 
with parameters a = v/2 and x = x 2 / 2 - 

The central gamma distribution functions satisfy the complementary re¬ 
lation 


P(a, x) + Q(a, x) — 1. (3) 

In hypothesis testing, and using the chi-square distribution, it is usual to 
consider the problem as Q(x 2 \ u ) = cr, a is the probability that a variable 
distributed according to the chi-square distribution with v degrees of freedom 
exceeds the value x 2 - Two kinds of problems are usually considered: a) 
computing the confidence level a given an experimentally determined chi- 
square value x 2 ; b) computing percentage points for certain values of a and 
for different degrees of freedom v (the standard tables which are common 
in statistics text-books usually include this information). The first problem 
involves the direct computation of the cumulative distribution function while 
the second one is an inversion problem. 

The module GammaCHI includes routines both for the direct compu¬ 
tation of the distribution functions (central gamma/chi-square) and their 
inversion. Regarding the direct computation, the routine first computes 
min{P(a, x), Q(a, x)}, and the other one by using (J3J) . Computing the Q(a, x) 
function simply as 1 — P(a,x) when P(a, x) is close to 1 can lead to serious 
cancellation problems, which are avoided in our algorithms. 

On the other hand, the inversion routine solves the equations 

P(a,x)=p, Q(a, x) = q, 0 < p, q < 1, (4) 

for a given value of a. 
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The algorithm for computing the gamma distribution function and its 
inversion is described in |6j]. The computation of the distribution functions 
is based on the use of Taylor expansions, continued fractions or uniform 
asymptotic expansions, depending on the parameter values. For the inver¬ 
sion, asymptotic expansions in combination with high-order Newton methods 
are used. 

The version included in the package GammaCHI for the inversion in¬ 
cludes some improvements in its performance for small values of p and q. 
For such small values, some of the elementary bounds for incomplete gamma 
functions given in [71] provide accurate enough approximations for estimating 
sufficiently accurate starting values for the Newton iterative process. 

2.2. Noncentral gamma distribution function 

The noncentral gamma cumulative distribution functions can be defined 
as 


Pn(x, y) = e + k > v), 




k=0 


where P(y,y) and Q(p,,y) are the central gamma cumulative distribution 
function and its complementary function, respectively; x is called the non¬ 
centrality parameter of the distribution functions. 

As in the case of the central distribution functions, the noncentral func¬ 
tions satisfy the relation 


Ppfay) + Qn(x,y) = 1. 


( 6 ) 


On the other hand, if XnW is a random variable with a noncentral chi- 
square distribution with n > 0 degrees of freedom and noncentrality param¬ 
eter A > 0, then the lower and upper tail probabilities for the distribution 
function of Xn( A) are given by 



k=0 




k=0 
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from which, and comparing with ([5]), it follows the simple relation between 
distribution functions: if the random variable Y has a noncentral gamma 
distribution with parameters a and noncentrality parameter A, then A" = 2 Y 
has a noncentral chi-square distribution with parameter n = 2a and with 
noncentrality parameter 2A. 

Algorithms for computing the functions P^(x,y) and QJ.x,y) for a large 
range of the parameters p (/r > 1), x, y are described in [8j. The methods 
include series expansions in terms of the incomplete gamma functions, recur¬ 
rence relations, continued fractions, asymptotic expansions, and numerical 
quadrature. In the module GammaCHI the lower /j range of computation 
of the algorithm for computing the noncentral gamma cumulative distribu¬ 
tion function is extended to include the interval | < p < 1. For extending the 
algorithms to this interval, we use a single step of the following three-term 
homogeneous recurrence relation (Eq.14 of jsj]) 

a , n , _ n _ [y h (2 y/xy) , , 

Ufi+i (1 + c^jy^ + c^y^-i 0, W — . (8) 

V x 4-1 (2 y/xy) 

Both Qn(x, y) and P^(x, y) satisfy (jSJ). In this expression, ratios of Bessel 
functions appear. These ratios are efficiently computed using a continued 
fraction representation. 

Regarding the inversion of the noncentral gamma/chi-square distribution 
functions, two different kinds of inversion problems can be considered: 

Problem 1: Find x from the equation 

Qn(x,y) = q, (9) 


with fixed y, /a and q. 

Problem 2: Find y from the equation 

P„(x, y) = p, 


( 10 ) 


with fixed x, /i and p. 

The inversion of Q^(x,y) with respect to x corresponds to the problem 
of inverting the distribution function with respect to the noncentrality pa¬ 
rameter given the upper tail probability. On the other hand, the inversion 
of P^(x,y) with respect to y with fixed x corresponds to the problem of 
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computing the p-quantiles of the distribution function. For noncentral chi- 
square distributions, the inversion with respect to y allows the computation 
of a variate with a given number of degrees of freedom and a given noncen- 
trality parameter. 

In the inversion with respect to x, it is important to note that not for all 
the possible values of the pair ( y 0 , qf) there is a solution of the equation 


Qv(x,yo)=qi- (11) 

In particular, the values of q x such that qi < q 0 , where q 0 = Q^(0, y 0 ) = 
Q[i{yo), do not make sense because Q^(x,y) is increasing as a function of x. 

Both in the inversion with respect to x and y, a combination of asymp¬ 
totic expansions and secant methods can be used, as discussed in |l[. A 
Fortran 90 implementation is now available in the module GammaCHI. As 
for the direct computation of function distribution values, in this implemen¬ 
tation the range of computation for the inversion of the gamma cumulative 
distribution function is extended to include the interval 1 < y < 1. The 
computation in this interval is made by using the double asymptotic prop¬ 
erty of the expansions derived in [lj: the formulas obtained in jll] can all be 
used for small values of y if £ = ‘lyjxy is large. More details can be found in 
the appendix. For small q this is suplemented with estimations for the case 
y = 1/2, which can be written in terms of error functions; we use the result 
of the inversion for the case /i = 1/2 to estimate a starting value when y is 
close to 1/2. It is easy to check that 

Qi/ 2 (x, y) = \ (erfc(v / x + y/y) + erfc (y/y - y/x)) 

and for this function, computing derivatives is simple and one can apply the 
fourth order method described in [9j]. 


2.3. Error function and its complementary function 

The error function erf (a:) jioj is defined by the integral 


erf(;r) = 


dt. 


7T Jo 


Its complement with respect to 1 is denoted by 


erfc (a:) = 1 — erf (a;) 



( 12 ) 


(13) 
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The complementary error function is very closely related to the normal 
distribution functions, which are defined by 

1 rx -t /»oo 

P(x) = ; _ / Q(x) — ._ / e~ t2 ^ 2 dt (14) 

V 2 W_oo \/ 2 W* 

with the property P(x) + Q(x) — 1. 

From the definition of the complementary error function it follows that 

P(x) = ^erfc(— x/V2), Q{x ) = -erfc(x/-\/ 2 ). (15) 

The algorithm for the error functions is based on the use of rational ap¬ 
proximations 0 and it includes the possibility of computing the scaled com¬ 
plementary error function e x erfc(x). This is important to avoid underflow 
problems in numerical computations of the complementary error function 
(the dominant exponential behavior of the function is e~ x ). In general, it is 
quite useful (when possible) to define scaled values of mathematical functions 
with the dominant exponential behavior factored out (Ti), §12.1.3]. 

Our module includes also the computation of the inverse of the comple¬ 
mentary error function y = inverfc x. Using the relation (|T5j) . this routine 
can be also useful in the generation of normal variates. 

For values of x close to 1, the following expansion is used 


inverfc x = t + t + + ..., 0 < x < 2 . 

6 30 630 7 


(16) 


with t = (1 — x). For these and more coefficients, see 13 


2-4- Gamma function 

Three standard definitions of the gamma function are usually considered: 


/*oo 

r(«) = / e~ l t z ~ l dt, Re z > 0 (Euler). (17) 

Jo 

Til Tl Z 

F(s) = lim ' — , 2 ^ 0 , 1, 2,.... (Euler). (18) 

n-too z[z + 1 ) ■ ■ ■ (z + n) 


1 

r® 


ze fZ jj (i + 

n=1 V 



e~ z,n 


(Weierstrass), 


(19) 


with 7 = 
proved in 


04)7721..., Euler’s constant. The equivalence of these definitions is 
14 . 










Important relations are 


r(z + l) =zT{z), (recursion), (20) 

and 


1 sin 7 rz 

T(1 + *)r(l — z) 7rz 


( 21 ) 


which is called the reflection formula. This formula is important for range 
reduction and it is easily proved by using (TTUT) . 

For computing the gamma function, we use recursion for small values of 
x and the relation given in (122|) for larger values. The special cases x = n 
and x = n + | are treated separately. 


2-4-1- Logarithm of the gamma function 

The following representation of logT(^) is very important for deriving 
expansions for large values of \z\: 

logT(^) = log z - ^ + \ log(2vr) + S(z). (22) 

For large \z\, S(z) gives a small correction with respect to the remaining 
terms of the right-hand side. Neglecting S(z) gives the well-known Stirling 
formula 


r(z) ~ z z e 


-(f) 


1/2 


z —> oo. 


(23) 


S(z) can be written in terms of an expansion involving Bernoulli numbers: 

N-l 


SM = E 


II 


2n+2 


—' (2 n + l)(2n + 2 )z 


2n+l 




(24) 


n =0 


The Bernoulli numbers B n are given by B n = B n ( 0). The first few are 


Bq — 1, B\ — —i?2 — B 3 — 0, B± — 


'311- 


More information on 


Bernoulli polynomials and Bernoulli numbers can be found in 15|, §24.2(i)]. 
A bound for in 


\Em\ < 


can be obtained as 

B 2 N +2 E(z) 

(2N + 1)(2AT + 2) Z 


-2N-1 


(25) 
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where 


K (z) = sup 

u>0 


2 i 2 

Z + U 


If |arg | < 7 t/4 , then K(z) = 1. For real positive z, En{z) is less in 
absolute value than the first term neglected in (1241) and it has the same sign. 
From (1251) it follows that E N (z ) = 0(z^ 2N ^ 1 ) for Re z —» +oo. Inserting the 
values of the first Bernoulli numbers we arrive at the representation (Stirling’s 
series): 


log r(z) = (z - ^) l°g £ - £ + \ log(27r) + 
1 


1 

TZz 


+ 


2 

1 

1260z 5 


+ 0(z~ 9 ). 


(26) 


360z 12601680^ 

The algorithm for the logarithm of the gamma function makes use of the 
representation (1221) . The computation of S(z) is based on the use of the sum 
(1221) . Chebyshev series and a minimax rational approximation. Coefficients 


of rational approximations are taken from the the tables in 16, § 6 . 6 ]. On 


some intervals we use Chebyshev expansions, with coefficients derived by 
using high precision Maple codes. 


2.5. The function T*(x) 

This function (the regulated gamma function) is defined by 


r*(*) 


r O) 

a/ 27 t/x x x e~ x ’ 


x > 0 . 


(27) 


When x is large, this function can be very important in algorithms where 
the function T(x) is involved because r*(x) = 1 + 0(l/x), as can be seen 
from its asymptotic expansion (Stirling series): 


r*(x)~l + ^a; 1 + -^x 2 + ..., x y oo. (28) 

The algorithm for computing T*(x), uses the relation T*(x) = exp(S'(a;)) 
for x > 3, where S(x) is given in (1241) . The computation of S(x) was described 
in the previous section. 
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2.6. The quotient of two gamma functions 

The quotient of two gamma functions appears frequently in applications. 
From a computational point of view, problems can arise when trying to 
compute directly the ratio of functions in the case when the arguments of 
both functions are large. For that reason, a key point in our algorithm for 
computing the quotient_of two gamma functions is the use of the following 
asymptotic expansion [l7j, p. 68]: 


r(x + a) 
r(x + b) 


w 


,a—b 


£(-i)"c. 


T(b-a + 2 n) 1 


n=0 


F(6- 


w 


.2 n ’ 


as x —» oo, 


with w = x +(a+ b — l)/2 and b > a. 

The first five coefficients C n appearing in (1291) are: 


Co = 1, 
r> — P 

“ 12 ’ 2 

C 2 = + 

r - p 1 p 1 p 

° 3 - gu72D + vrm + id3bs> 

r - p 1 mi c 1 p 1 p 

° 4 - 4838400 + iUi 87091200 + 414720 + 497664 ’ 
with p — (a — b + l)/2. 


(29) 


(30) 


3. Overview of the software structure 

The Fortran 90 package includes the main module GammaCHI, which 
includes as public routines the following functions and routines: cdfgamC, 
invcdfgamC, cdfgamNC, invcdfgamNC, errorfunction, inverfc, gamma, 
loggam, gamstar and quotgamm. 

In the module GammaCHI, the auxiliary module Someconstants is 
used. This is a module for the computation of the main constants used in 
the different routines. 


4. Description of the individual software components 
1. cdfgamC: 

The calling sequence of this routine is 
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cdfgamC(ichi,a,x,p,q,ierr) 


where the input data are: ichi (flag), a and x (arguments of the distri¬ 
bution function), ichi is a flag for the choice of a gamma or chi-square 
cumulative distribution function: when ichi = 1, the arguments are of 
a gamma distribution function; when ichi = 2, the arguments corre¬ 
spond to a central chi-square distribution function. 

The outputs of the function are error flag ierr, the function distribution 
value p and its complementary function q. The possible values of the 
error flag are: ierr = 0, successful computation; ierr = 1, computation 
failed due to overflow/underflow; ierr = 2, arguments out of range. 

2. invcdfgamC: 

The calling sequence of this routine is 
invcdfgajnC(ichi ,a,p,q,xr, ierr) 

where the input data are: ichi (flag), a, p and q. ichi is a flag for the 
choice of a gamma or chi-square cumulative distribution function: when 
ichi — 1, a gamma distribution function is inverted; when ichi = 2, a 
central chi-square distribution function is inverted. 

The outputs of the function are xr (the solution of the equations 
P(a,xr) = p and Q(a,xr) = q) and the error flag ierr. The pos¬ 
sible values of the error flag are: ierr = 0, computation successful; 
ierr = 1, overflow problems in one or more steps of the computation; 
ierr = 2, the number of iterations in the Newton method reached the 
upper limit N = 35; ierr = 3, any of the arguments of the function is 
out of range. 

3. cdfgamNC: 

The calling sequence of this routine is 
cdfgamNC(ichi,mu,x,y,p,q,ierr) 

where the input data are: ichi (flag), mu, x and y (arguments of the 
distribution function). The input ichi input is a flag for the choice 
of a gamma or chi-square noncentral cumulative distribution function: 
when ichi = 1, the arguments are of a gamma distribution function; 
when ichi = 2, the arguments correspond to a noncentral chi-square 
distribution function. 


12 



The outputs of the function are: the function distribution value p and 
its complementary function q and an error flag ierr. The possible 
values of the error flag are: ierr = 0, successful computation; ierr = 1, 
computation failed due to overflow/underflow; ierr = 2, arguments out 
of range. 

4. invcdfgamNC: 

The calling sequence of this routine is 

invcdfgamNC(ichi,icho,mu,p,q,yx,xy,ierr) 


where the input data are: ichi, icho, mu, p, q and yx. ichi is a flag for 
the choice of a noncentral gamma or chi-square cumulative distribution 
function: when ichi = 1, the arguments are of a gamma distribution 
function; when ichi = 2, the arguments correspond to a central chi- 
square distribution function. The input icho is a flag for the choice 
of the inversion process: if icho — 1 , x is computed in the equations 
Qmu( x ,yx) = q, P m u( x ,y x ) = V • If icho — 2, y is computed in the 
equations Qmu(u x ,y) = Q, Pmu(y x ,y) — P', Q and P are the values of 
the cumulative distribution functions; mu and yx are arguments of the 
distribution functions. 

The outputs of the routine are the value obtained in the inversion 
xy and an error flag ierr. The possible values of the error flags are: 
ierr = 0, computation successful; ierr = 1, in the inversion process 
with respect to the non-centrality parameter {x- variable), the input 
values (y, q ) do not make sense (i.e. the inequality q < Q mu ( 0, y) holds); 
ierr = 2, at least one of the gamma distribution function values needed 
in the inversion process cannot be correctly computed; ierr = 3, the 
number of iterations in the secant method reached the upper limit; 
ierr = 4, any of the input arguments is out of range. 

5 . err or function: 

The calling sequence of this function is 
errorfunction(x,erfc,expo) 


where the input data are: x (argument of the function), erfc and expo 
(logical variables). 

The logical variables erfc and expo mean the following: 
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When erfc=.true. and expo=.false., the function computes the com¬ 
plementary error function erfc(x). 

When erfc=.true., expo=.true. and x > 0, the function computes the 
scaled complementary error function e x erf'c(x). 

When erfc=.false, and expo=.false., the function computes the error 
function erf(x). 

6. inverfc: 

The calling sequence of this function is 
inverfc(x) 

The input argument is x, a positive real number. This function com¬ 
putes the inverse of the complementary error function. 

7. gamma: 

The calling sequence of this function is 
gamma(x) 

This function computes the Euler gamma function T(x), for x real. 

8. loggam: 

The calling sequence of this function is 
loggam(x) 

This function computes log(T(x)), for x > 0. 

9. gamstar: 

The calling sequence of this function is 
gamstar(x) 

This function computes for positve x the regulated gamma function 

r*(x) = . r (~ r ). _ . 

\/2n/x x x e x 

10. quotgamm: 

The calling sequence of this function is 
quotgamm(x,y) 


This function computes 


EM 

rfe) 


for x, y real values. 
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a 

X 

Q(a, x) 

gammainc(x, a,' upper 1 ) 

IQ-250 

6.3 10 -15 

3.212101109661167 10~ 249 

3.21210110966116710~ 249 

10-25° 

7.110~ 7 

1.3580785912009393 10~ 249 

1.35807859120093 10" 249 

10-25° 

0.01 

4.0379295765381135 10~ 25 ° 

4.037929576538114 10~ 25 ° 

IQ- 14 

6.3 10 -15 

3.21210110966065110~ 13 

3.21210110966065210" 13 

IQ- 14 

7.110 -7 

1.358078591200848 10~ 13 

1.35807859120084810" 13 

IQ- 14 

0.01 

4.0379295765380405 10" 14 

4.037929576538039 10~ 14 


Table 1. Test of the routine cdfgamC for small values of a and x: comparison of 
the values obtained for Q(a,x ) using our routine cdfgamC and the Matlab function 
gammainc included in the 2013a release. It should be noted that our value for x = 0.01, 
a = 10“ 14 agrees in 16 digits with the value obtained with the Mathematica function 
GammaRegularized[a,x]: 4.03792957653804043... 10 -14 . 


5. Testing the algorithms 


For testing the algorithm for the computation of the central gamma dis¬ 


tribution functions, we use the recurrence relations (see 18, §8.8]) 


P(a + 1, x) = P(a, x ) — D(a, x ), 


Q(a + 1, x) = Q(a, x) + D(a, x), (31) 


where 


D(a, x) 


r(a + 1 ) 


(32) 


The recurrence relations (l3Tj) are tested for a large number of random 
points in the input parameter domain (a, x), x < 0, a > 10~ 300 . The aimed 
accuracy is close to 10 -13 in most cases. For small values of a, the accuracy 
of the values computed by our routine cdfgamC is even higher, near 10 15 
as can be seen in Table 1. In this table, a comparison of a few values against 
those obtained using the Matlab function gammainc included in the release 
2013a, is shown. It should be noted that versions of gammainc included in 
some previous Matlab releases showed some loss of accuracy for these same 
values. 

For the inversion algorithm, testing is made by checking that the compo¬ 
sition of the functions with their inverse is the identity. Some tests for the 
inversion routine invcdfgamC are shown in Table 2. This table shows the 
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a 


p 

0.05 

1 

10 

100 

1000 

0.0001 

4 10 -lti 

4 . 110 " ie 

6.5 10“ 15 

8 . 110 -lti 

5.3 10" 15 

0.01 

2 10 -16 

5.2 10~ 16 

2.9 10" 15 

1.2 HT 15 

1.9 10" 15 

0.1 

2.8 10" 16 

1.4 10 -16 

0 

1.5 10~ 15 

8.3 10" 16 

0.3 

0 

1.8 10 -16 

3 . 710 -16 

7.4 10~ 16 

7.4 10~ 16 

0.5 

io - 16 

0 

0 

1 . 110~ 16 

2.2 10~ 16 

0.7 

0 

1.6 10 -16 

0 

1.6 10~ 16 

7.9 10" 16 

0.9 

0 

0 

0 

0 

3 . 710~ 16 

0.9999 

0 

0 

0 

0 

0 


Table 2 . Relative errors \P(a,x r ) — p\/p for several values of p and a. 


relative error for the inversion problem of the central gamma distribution 
function given in Eq. (JU) for different values of p and a. 

Further tests for the central gamma distribution functions can be found 
in |6]. 

Testing of the noncentral gamma distributions can be made by checking 
the deviations from 1 of the relation 

(x - p)Q^ +1 (x, y) + (y + p)Q^(x, y) _ ^ 

xQi*+2{x,y) + yQ»-i(x,y) 

Equation (1331) is used when y > x + p, and the same expression but for 
P, (x, y ) when y < x + y. The aimed accuracy for the computation of the 
noncentral gamma cumulative distribution functions is close to 10~ n in most 
cases in the parameter domain (x, y, /i) with 0 < x < 10000, 0 < y < 10000 
and 0.5 < p, < 10000. 

We have compared our routines for the noncentral distributions against 
the package ncg in R. The description of the algorithms included in this pack¬ 
age is given in [19| . Algorithms (but not associated software) for the compu¬ 
tation of the noncentral gamma distributions are also described in 0. The 
package ncg includes functions for the computation of the noncentral gamma 
cumulative distribution function (function pgammanc), the computation of 
the noncentrality parameter (function deltagammanc) and the quantiles of 
the noncentral gamma distribution function (function qgammanc). In our 
notation, this corresponds to the computation of the function P A1 (x,y), the 
inversion of P M (x, y) with respect to x and the the inversion of P At (x, y) with 
respect to y, respectively. Therefore, there is a first difference between our 
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/i, x j y 

pgammanc(t/, /i, 2x) 


5, 150, 30 

1.132428676150... 10~ 23 

GST: 1.215915354045...10 23 

MATH: 1.215915354045. ,10” 23 

1, 75, 0.5 

1.215915354045...10" 23 

GST: 3.287840255874...10 30 

MATH: 3.287840255874...10- 30 

2, 100, 2 

1.976996751829... 10 -38 

GST: 1.557081489535. ..10 35 

MATH: 1.557081489535... 10 -35 

10, 100, 1 

5.486504449527...10- 57 

GST: 5.152185145235...10 48 

MATH: 5.152185145235...10 48 


Table 3. Test of the routine cdfgamNC: comparison of the values obtained for P^x^y) 
using the R function pgammanc, our routine cdfgamNC (values GST) and the Matli- 
ematica function 7V[MarcumQ[/u, y/2x, 0, \/2y]] (values MATH) for some values of /i, x 
and y. 


routines and those included in the package ncg: our algorithms for the direct 
computation allow to obtain both P^x^y) and its complementary function, 
the function Q^Xjy). Also, the inversion algorithms compute for the mini¬ 
mum of P /i (a;,t/) and Q^(x, y), which is very convenient for computations in 
the queues of the distributions. 

In the computation of P^(x, y), we have checked that the function pgam¬ 
manc has some problems when computing small values of the cumulative 
distribution function. This is illustrated in Table 3, where computations by 
pgammanc, our routine cdfgamNC and the Mathcmatica function Mar- 
cumQ for some values of /i, x and y are shown. 

Regarding the comparison of the inversion routine with respect to x (the 
noncentrality parameter), the computations using the function deltagam- 
manc included in the package ncg seem to be very slow for certain pa¬ 
rameter regions (small p, large y, small p). For example, for computing x 
in Pi.g(a;,288) = 0.00001 (2x when calling the function deltagammanc) 
was about 40 s when using deltagammanc and 10” 5 s when using invcd- 
fgamNC. 

Another test for the inversion of the noncentral gamma distribution func- 
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Inv. wrt x 

Inv. wrt y 

q 

y = io 

y = 100 

y = 1000 

x = 10 

x = 100 

x = 1000 

0.001 

4 10 -13 

4 10~ ie 

1.5 10~ 15 

6.5 10 -lti 

1.5 10" 15 

1.5 10 -15 

0.1 

2 10~ 15 

810" 14 

4 10 -14 

2.8 10 -16 

8.110" 14 

2.5 10 -13 

0.3 

1IO" 13 

0 

3 10 -15 

9.3 10 -16 

1.5 IO" 14 

7.3 10 -14 

0.5 

4 10 -15 

0 

0 

2.5 10~ 14 

0 

2 10 -15 

0.7 

5 10~ 15 

8 10 -16 

4.9 10 -15 

5 10 -15 

3.2 10" 14 


0.999 

2 10 -16 

0 

0 

0 

0 

0 


Table 4. Relative errors in the inversion of the noncentral gamma distribution function 
with respect to (wrt) x and with respect to y. The value of y has been fixed to 0.5 


tion, is shown in Table 4. In this table, the relative error for the inversion 
problem with respect to x and y of the noncentral gamma distribution func¬ 
tion given in Eq. (0J for different values of (q, x ) and (q, y ) (/i fixed to 0.5) 
is shown. In a second test, 10' random points have been generated in the 
domain of input parameters fi G [0.5, 10000], y G [0, 10000] (x G [0, 10000]) 
and q G [0, 1] for testing the inversion with respect to x (y). Fixing the accu¬ 
racy parameter of the secant method to 10~ 12 , the accuracy of the computed 
values obtained in the inversion was always better than 10~ n . 

The accuracy of the routines for the computation of the error function, the 
complementary error function, the gamma function and its related functions 
have been tested by comparing the values obtained against Maple (with a 
large number of digits). In all cases, a full agreement in at least 13 — 14 
significant digits was found. 

Finally, Table 5 shows relative errors in the inversion process of the com¬ 
plementary error function y = erfc(x) for different values of y. 

6. Test run description 

The Fortran 90 test program testgamCHI.f90 includes several tests for 
the computation of function values and their comparison with the corre¬ 
sponding pre-computed results. Also, the inversion routines are tested by 
checking that the composition of the function with its inverse is the identity. 
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y 

X 

||erfc(a;)/t/ — 1|| 

1.9 

-1.1630871536766738 

0 

1.0 

-8.571414913863924 10" 17 

0 

nr 2 

1.163087153676674 

1.4 10 -16 

10- 3 

2.3267537655135246 

4.3 10“ 16 

10- 4 

2.7510639057120607 

2.710“ 16 

10- 5 

3.123413274340875 

1.2 10“ 16 

10“ 6 

3.4589107372795 

8.4 10" 16 

10- 7 

3.766562581570838 

6.6 10 -15 

10“ 8 

4.052237243871389 

2.110 -15 

10- 9 

4.320005384913445 

8.3 10 -16 

10- 10 

4.572824967389486 

4.8 10- 15 

10- 11 

4.812924067365823 

4.8 10 -16 

10- 12 

5.042029745639059 

0 


Table 5. Solutions x(y) of the equation y = erfc(a;). 


7. Appendix: The asymptotic inversion method 

The cumulative distribution functions considered in this paper can be 


written in one of the standard forms (see also 21 



1^1 e-*°< 


s? 




(34) 


where a > 0, J) 6 K, and / is analytic and real on M with /(0) = 1. The 
special case / = 1 gives the normal distributions (see (fTTD ) 


P(vVa) = J- ^ J e ^ d( = |erfc (^rjy/aj^j , 

QivVa) = f e ~^ <2 dC > = 2 erfc (+^^/2) • 

We explain how the incomplete gamma function P(a,x ) defined in 
can be written in this standard form. Let A = - and t = ar. Then 


(35) 


T*(a)P(a, x ) = 



a 

27T 


3 —a(r—lnr—1) 


dr 


T 


(36) 
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( 37 ) 


where T*(a) is defined in (1271) . The transformation 

r - lnr - 1 = |C 2 , sign(r - 1) = sign(C), 
gives the standard form 

r(a)P(a,x) = f /(C) = (38) 

where r; is defined by 

A — In A — 1 = 17? 2 , sign(A — 1) = sign(r/). (39) 


By using an integration by parts procedure in we can obtain an 
asymptotic representation of F a (rf) which is valid for a —» oo, uniformly for 
all rj £ M. We write f(rj) = [fin) — /(0)] + /(0), where /(0) = 1, and use 
(l35]h Then we obtain by repeating integration by parts 


_ e ~2 ar l 

Fair]) = |erfc(-r 7 V / o 72 )^a(oo) + — S a {rj), 

\/27ra 

,_ e ~k ar i 2 

Gain) = §erfc(?/v/a72)G a (—oo) - — S a (ri), 

V 2vra 

where F a (oo) = G a (—oo) and 


(40) 


A 


F ai OO) ~ -?> A) = 1, AW ~ 


a(//) 


a —> oo, 


(41) 


n=0 


n=0 


uniformly with respect to p £ M. The coefficients follow from the following 
recursive scheme. Let fo(n) = fin)- Then, for n — 0,1, 2,..define 


fn+i in) 


d fnjy) - /»(0) 

dn n 


(42) 


and we have 

A, = /„(0), C„fo) = (43) 

V 

To describe the inversion problem, we assume that p £ (0,1) and that 
a is a large positive parameter. Then we are interested in the value n tliat 
solves the equation 

F a in) — F a {oo)p, (44) 
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where F a (rj) has the form given in (TUI . 

First we define a number that solves the reduced equation 

|erfc (-rj 0 y/a/^Zj = p. (45) 

Then for the desired value 77 we assume the expansion 

h ~ ho + — + “T + ^| + ..., a-> oo, (46) 

a a z cr 

and try to find the coefficients r/i, r/ 2 , r/ 3 ,.... To obtain the 77 ^ we can sub¬ 
stitute the expansion for 77 into the asymptotic expansion of F a {rj) and use 
formal power series manipulations. 

However, here we use the method as explained in our earlier mentioned 
publications. This method runs as follows. From (l34]h (I44H . and (1451) we 
obtain _ _ 

4 [± e -wo h = A r* m. c -w 

drj 0 V 2vr ’ dr] y 27rF a (oo) 

from which we obtain, upon dividing, 

f(ri)F = Fci{ao)e H^\ (48) 

Substituting (1461) and using F a (oo ) = 1 + C?(l/a), we obtain, after hrst-order 
perturbation analysis for large a, 

f(r ?o) = e. VoVl => rji = — 111 /( 770 ). (49) 

Vo 



Because / is analytic at the origin with /(0) = 1, ifo is well defined as r] 0 —* 0. 

For higher-order terms tj 3 , 7 > 2, we need in (|48|) more coefficients in the 
asymptotic expansion of F a (oo ) (see (14B ) and we have to expand 

/(h) = /(ho) + (h ~ ho)/'(ho) + |(h - ho) 2 /"(ho) + • • • • (50) 


Then, the next coefficients are given by 

h2 = - (/ (2^i + / - 2p() - 2771 /) / (2p 0 /), 

h3 = - (8/hih2 + 4A!/?7/ + frjf + 4/77/770772 + 4/77^77! - 8/77^+ (51) 

8A1/770772 + 8H 2 / - 8/77177^ - 8/772 - 4/V) / (8770/), 
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where /, f and f" are evaluated at p 0 . 

For small values of (that is, when p ~ |), we need expansions. We 
have 

Vi = ai + \ (2 a 2 - a?) 77o + | (3a 3 - 3aia 2 + 2a?) r/„ + C> ( 7 $) , 

772 = — |n‘i + 2a 3 + | ( —12a 2 cii + 5a^ + 24a 3 ) r/o + (9 ( 770 ) , (52) 

V 3 = pg (4af - 5a 2 ai - 15a 3 cii + 120a 5 ) + C> ( 7 / 0 ), 


where are the coefficients in the expansion /(£) = ^^ak( k - The afc are 

k =0 

related to the A k used in (T5T]l : 

(53) 


A n =(±) 2 n a 2n , 77 = 0,1,2,,.... 

\ / n 
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